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Abstract 


In  this  paper,  several  numerical  aspects  of  an  existing  model  for  fully  nonlinear  waves  are 
improved  and  validated  to  study  wave  breaking  due  to  shoaling  over  a  gentle  plane  slope  and 
wave  breaking  induced  by  a  moving  lateral  boundary. 

The  model  is  based  on  fully  nonlinear  potential  flow  theory  and  combines  a  higher-order 
Boundary  Element  Method  (BEM)  for  solving  Laplace’s  equation  at  a  given  time  and  Lagrangian 
Taylor  expansions  for  the  time  updating  of  the  free  surface  position  and  potential.  An  improved 
numerical  treatment  of  the  boundary  conditions  at  the  intersection  between  moving  lateral 
boundaries  and  the  free  surface  (comer)  is  implemented  and  tested  in  the  model,  and  the  free 
surface  interpolation  method  is  also  improved  to  better  model  highly  curved  regions  of  the  free 
surface  that  occur  in  breaking  waves.  Finally,  a  node  regridding  technique  is  introduced  to 
improve  the  resolution  of  the  solution  close  to  moving  boundaries  and  in  breaker  jets. 

Examples  are  presented  for  solitary  wave  propagation,  shoaling,  and  breaking  over  a  1:35  slope 
and  for  wave  breaking  induced  by  a  moving  vertical  boundary.  Using  the  new  methods,  both 
resolution  and  extent  of  computations  are  significantly  improved  compared  to  the  earlier  model, 
for  similar  computational  efforts.  In  all  cases  computations  can  be  carried  out  up  to  impact  of 
the  breaker  jets  on  the  free  surface. 
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1  Introduction 


Over  the  past  fifteen  years,  efficient  models  were  developed  for  the  solution  of  two-dimensional 
fully  nonlinear  potential  flows  with  a  free  surface  (FNPF)  (e.g.,  Longuet-Higgins  and  Cokelet 
(1976)  (LHC);  Vinje  and  Brevig  (1981);  New,  Mclver  and  Peregrine  (1985);  Dold  and  Peregrine 
(1986)  (DP);  Grilli,  Skourup,  and  Svendsen  (1989)  (GSS)).  A  few  attempts  have  also  been 
reported  in  solving  FNPF  problems  in  three  dimensions  (Romate  (1989)  (ROM);  Xii  and  Yue 
(1992);  Broeze  (1993))  (see  also  the  reviews  in  (GSS;  ROM).  In  most  successful  approaches, 
Laplace’s  equation  is  solved  at  a  given  time  based  on  a  Boundary  Integral  Equation  (BIE)  or 
on  a  Boundary  Element  Method  (BEM),  and  free  surface  geometry  and  potential  are  updated  in 
time  based  on  a  mixed  Eulerian-Lagrangian  description  of  the  free  surface. 

In  the  earlier  models,  waves  were  assumed  to  be  periodic  in  space  and  the  problem  was  then 
solved  in  a  transformed  space  and/or  using  a  complex  variable  formulation  (e.g.,  LHC;  DP). 
To  apply  the  method  to  arbitrary  waves,  complex  bottom  topography,  or  moving  boundaries, 
other  models  were  proposed  which  solved  the  equations  directly  in  the  physical  space  (Lin, 
New  and  Newman  (1984)  (LNN);  GSS;  Cointe  (1990)).  [Note  also  that,  for  cases  in  which 
the  unsteady  flow  vanishes  in  both  horizontal  directions,  Cooker  (1990)  extended  the  space 
periodic  model  by  DP  to  irregular  bottom  topography.]  These  models,  however,  sometimes 
had  problems  when  simulating  the  generation  of  waves  by  wavemakers  (e.g.,  LNN).  Among 
these  early  models,  the  model  by  GSS  and  Grilli  and  Svendsen  (1990)  (GS),  was  the  object  of 
considerable  developments  and  improvements  which  made  it  particularly  versatile  and  accurate 
for  solving  problems  of  wave  breaking  induced  by  fixed  or  moving  boundaries.  In  this  model, 
Laplace’s  equation  is  solved  by  a  higher-order  BEM  based  on  2nd  Green’s  identity  (a  method 
also  applicable  to  three-dimensional  problems  unlike  earlier  approaches  based  on  complex 
variables),  using  isoparametric  elements  for  the  spatial  discretization  (2nd  to  4th-order),  and 
time  integration  is  explicit  and  based,  as  in  DP,  on  (2nd-order)  truncated  Taylor  series  expansions. 
More  details  about  this  numerical  model  can  also  be  found  in  Grilli  (1993). 

When  working  in  the  physical  space,  incident  waves  must  be  generated  and  reflected  or 
absorbed  in  the  model,  which  poses  much  more  difficult  problems  than  with  space-periodic 
models.  The  modeling  of  wave  absorption  using  a  “sponge  layer”  is  discussed,  e.g.,  in  Cointe 
(1990)  and  Grilli  and  Horrillo  (1995)  (see  Otta,  Svendsen  and  Grilli  (1992)  (OSG)  for  a  literature 
review).  The  modeling  of  the  intersection — or  comer — between  the  free  surface  and  a  fixed 
or  a  moving  boundary  has  been  the  object  of  a  substantial  amount  of  work  in  the  literature,  in 
the  context  of  wave  generation  by  a  wavemaker  starting  from  a  state  of  rest  (Roberts  (1987); 
Cointe  (1989);  GS  (also  for  a  literature  survey  to  date);  Joo,  Schultz  and  Messiter  (1990);  OSG; 
Svendsen,  Otta  and  Grilli  (1992)).  It  was  thus  found  that  the  sudden  jump  to  a  finite  value  of  the 
wavemaker  initial  velocity  and  acceleration  may  create  an  initial  (mathematical)  singularity  of 
the  flow  near  the  wavemaker.  This  singularity  results  from  inconsistencies  between  boundary 
conditions  on  the  free  surface  and  on  the  wavemaker  and,  if  not  properly  accounted  for  in  the 
numerical  model,  invariably  worsens  through  time  marching,  eventually  leading  to  instability  of 
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the  numerical  solution.  Such  a  problem  may  be  addressed  two  ways  in  a  numerical  model :  (i) 
either  by  explicitly  introducing  the  singular  behavior  of  the  flow  in  the  numerical  solution  (as, 
e.g.,  discussed  in  Cointe  (1988));  (ii)  or,  more  simply,  by  specifying  initial  boundary  conditions 
in  such  a  way  that  both  inconsistencies  and  singularities  do  not  occur  (OSG).  The  latter  method 
is  further  discussed  in  the  paper  and  will  be  used  in  the  present  applications. 

Besides  the  possible  occurrence  of  an  initial  (mathematical)  singularity  of  the  solution,  it 
can  be  seen  that  additional  problems — most  of  these  being  generic  to  BEM/BIE  methods — can 
occur  at  corners  as  a  result  of  improper  numerical  representation  of  the  solution.  In  fact,  without 
special  treatments  of  comers,  BEM/BIE  solutions  usually  exhibit  large  errors  at  comers,  even 
when  the  actual  mathematical  solution  is  non-singular  (e.g.,  Schultz  and  Hong  (1989)).  In  the 
context  of  wave  generation  by  wavemakers,  a  poor  numerical  solution  at  corners  may  thus  have 
all  the  appearances  of  singular  behavior  and,  through  time  marching,  lead  to  instability  of  the 
numerical  solution.  GSS  and  GS  analyzed  the  causes  for  such  (numerical)  singularity  problems 
and  proposed  specific  numerical  treatments  to  eliminate  them  :  (i)  the  representation  of  corners 
by  double-nodes  and  the  use  of  continuity  conditions  for  the  potential;  (ii)  the  specification  of 
fluid  velocity  uniqueness  at  comers  through  using  compatibility  conditions  between  the  normal 
and  tangential  derivatives  of  the  potential;  (iii)  the  improved  accuracy  of  (non-singular  but 
potentially  quasi-singular)  numerical  integrations  at  comers,  through  using  adaptive  numerical 
integration  techniques.  Using  such  treatments,  GS  calculated  the  solution  of  various  mixed 
boundary  value  problems,  in  computational  domains  with  sharp  corners,  with  almost  arbitrarily 
small  errors  and  GS;  Svendsen  and  Grilli  (1990);  Grilli,  Losada  and  Martin  (1994);  Grilli, 
Subramanya  and  Svendsen  (1994),  successfully  computed  wave  generation  by  a  wavemaker 
and  wave  propagation  and  runup  over  slopes  and  over  other  bottom  discontinuities  or  structures. 
Finally,  based  on  the  same  methods  as  in  GS,  OSG  developed  extended  corner  compatibility 
conditions  (numerical  treatment  (ii))  which  provided  further  improvements  of  the  corner  solution 
in  the  case  of  wavemaker  motion.  These  extended  conditions  are  detailed  in  the  present  paper 
and  will  be  tested  in  the  applications  versus  the  earlier  approaches  by  GSS  and  GS. 

In  the  model  by  GSS  and  GS,  as  in  most  earlier  models  based  on  a  mixed  Eulerian- 
Lagrangian  description,  free  surface  discretization  nodes  represent  water  particles  whose  posi¬ 
tion  is  calculated  as  a  function  of  time.  Doing  so,  multiple-valued  free  surface  elevations  can 
be  represented  and  overturning  waves  can  in  principle  be  calculated  up  to  the  time  the  tip  of  the 
breaker  jet  impacts  the  free  surface  (thus  violating  continuity  equation  and  leading  to  numerical 
instability  of  the  computations).  This  is  illustrated  in  the  applications  in  this  paper.  It  is  worth 
pointing  out  that  laboratory  experiments  and  computations  by  Dommermuth,  Yue,  Lin,  Rapp, 
Chan  and  Melville  (1988),  for  deep  water  breakers,  and  by  Grilli,  Subramanya,  Svendsen  and 
Veeramony  (1994)  (GSSV),  for  shallow  water  breakers,  both  confirmed  the  physical  validity  of 
the  FNPF  results  up  to  the  point  of  impact  (see  also  results  by  Grilli,  Svendsen  and  Subramanya 
(1994)). 

In  using  the  model  by  GSS  to  compute  overturning  waves,  GS  reported  numerical  inaccu¬ 
racies  of  the  solution  in  highly  curved  regions  of  the  free  surface,  like  breaker  jets,  before  wave 
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overturning  could  be  fully  calculated.  They  showed  that  these  inaccuracies  were  due  to  small 
discontinuities  in  the  free  surface  slope  occurring  between  successive  isoparametric  elements 
and  could  be  eliminated  with  a  representation  of  the  free  surface  geometry  ensuring  continu¬ 
ity  of  the  slope.  This  representation  used  linear  shape  functions  to  discretize  field  variables 
and  cubic  splines  to  represents  the  geometry,  and  was  referred  to  as  “quasi-spline  boundary 
elements”  (QS).  Despite  the  linear  shape  functions,  QS  elements  proved  quite  accurate  for  a 
wide  range  of  applications  of  the  model  (e.g.,  Grilli  (1993);  Grilli,  Losada  and  Martin  (1994)). 
A  more  detailed  analysis  of  this  problem  by  OSG  also  showed  that,  when  using  higher-order 
isoparametric  elements  on  the  free  surface — and  not  linear  elements  as  in  most  earlier  models 
(e.g.,  LHC;  VB) — due  to  the  shape  function  representation  of  the  geometry,  time  updating  of 
intermediate  nodes  in  the  elements  is  treated  slightly  differently  than  time  updating  of  extremity 
nodes.  Such  differences  in  the  updating  may  make  small  discontinuities  in  inter-element  slope 
grow  with  time,  particularly  in  regions  of  high  curvature,  and  trigger  larger  scale  instabilities  of 
the  free  surface.  In  fact,  using  simpler  piecewise  linear  elements  on  the  free  surface,  which  do 
not  have  intermediate  nodes,  GSS  were  able  to  compute  large  scale  plunging  breakers  identical 
to  those  reported  in  (LHC;  VB;  NMP)  without  facing  instability  problems  discussed  above. 
Overall,  however,  their  solution  was  less  accurate  than  when  using  higher-order  elements  that 
do  not  have  discontinuity  of  the  slope  (see  below).  This  is  because  a  higher-order  representation 
of  the  free  surface  improves  both  the  accuracy  and  the  convergence  of  the  BEM  solution  of 
Laplace’s  equation  (GSS;  OSG). 

OSG  thus  introduced  a  higher-order  representation  of  the  free  surface  for  both  the  geometry 
and  the  field  variables  that  ensures  continuity  of  the  free  surface  slope.  To  do  so,  they  defined  the 
interpolation  within  each  free  surface  segment  between  two  nodes  using  the  middle  interval  of 
a  four-node  cubic  isoparametric  element  (thereafter  referred  to  as  “mid-interval  interpolation” 
or  Mil)  \  In  the  computations,  the  Mil  element  is  slid  along  the  free  surface,  one  node  at  a 
time,  thus  defining  a  full  cubic  interpolation  (unlike  with  the  QS  elements)  with  a  continuous 
slope.  Using  the  Mil  method,  OSG  showed  improvements  of  several  orders  of  magnitudes  in 
computational  accuracy,  as  compared  to  the  QS  method,  without  significant  change  in  computing 
time.  Although  the  Mil  method  proved  very  accurate  for  the  solution  of  problems  with  fixed 
boundaries,  it  was  found  less  accurate  or  even  unstable  for  problems  with  strong  interactions 
between  a  moving  lateral  boundary  and  the  free  surface.  In  this  paper,  it  is  shown  that  results 
can  be  improved  for  such  cases  using  a  modified  Mil  method  in  which  a  cubic  spline  is  used  for 
the  interpolation  of  the  geometry  and  the  same  sliding  cubic  polynomial  as  in  the  Mil  method 
is  used  for  the  interpolation  of  the  field  variables.  This  new  interpolation  method,  thereafter 
referred  to  as  “mixed  cubic  interpolation”  or  MCI,  is  detailed  in  the  paper  and  computational 
examples  are  presented  for  wave  breaking  induced  by  fixed  or  moving  boundaries. 

Finally,  when  generating  waves  by  a  moving  boundary  (wavemaker),  it  is  observed  that  free 
surface  nodes  gradually  drift  away  from  the  moving  boundary  and  concentrate  in  convergence 
regions  of  the  flow  like  breaker  jets  (e.g.,  GSS).  This  node  drift  is  due  to  the  mean  nonlinear  mass 

3  This,  in  fact,  followed  the  same  principle  as  used  in  DP  for  space-periodic  problems. 
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flux  of  the  flow — the  so-called  Stokes  drift — and  is  predicted  by  nonlinear  wave  theories.  In 
the  present  case,  with  time,  node  drift  leads  to  a  poor  resolution  of  the  discretization  close  to  the 
wavemaker,  a  region  of  the  flow  with  large  variations  of  field  variables  requiring  good  resolution 
of  the  discretization,  and  to  a  high  concentration  of  nodes  in  breaker  regions,  which  may  create 
quasi-singular  situations.  Both  of  these  node  drift  effects  result  in  increased  computational 
errors,  which  are  discussed  in  the  paper,  and  an  adaptive  regridding  technique  is  introduced 
to  redistributes  nodes  along  a  numerically  calculated  free  surface  and  improve  the  accuracy  of 
computations  in  such  situations  4. 


2  The  mathematical  and  numerical  model 

2.1  Governing  equations  and  boundary  conditions 

Equations  for  the  two-dimensional  potential  model  by  (GSS;  GS)  are  briefly  listed  in  the 
following.  The  velocity  potential  <j>(x,  t )  is  used  to  represent  inviscid  irrotational  2D  flows  in 
the  vertical  plane  (x,z)  and  the  velocity  is  defined  by,  u  =  V</>  =  ( u,w )  (Fig.  1).  Continuity 
equation  in  the  fluid  domain  f l(t)  with  boundary  T(t)  is  a  Laplace’s  equation  for  the  potential, 

V2<j>  =  0  infi(i)  (1) 


Using  the  free  space  Green’s  function,  G(x,  x{)  =  — (l/27r)log  |  x  —  xt  |,  and  Green’s 
second  identity,  equation  (1)  transforms  into  the  Boundary  Integral  Equation  (BIE), 

a{xi)(j){xi)  =  xi )  “  </>(x)dG^Xl^}  dV{x)  (2) 

in  which  x  =  (x,  z)  and  xi  =  (xi,  z{)  are  position  vectors  for  points  on  the  boundary,  n  is  the 
unit  outward  normal  vector,  and  a{x{)  is  a  geometric  coefficient  function  of  the  exterior  angle 
of  the  boundary  at  aq. 


On  the  free  surface  T f(t),  <f>  satisfies  the  nonlinear  kinematic  and  dynamic  boundary 
conditions, 

^L  =  u  =  \;(j)  on  Tf(t)  (3) 

^7  =  -S*  +  \  ^  on  T,(t)  (4) 

U  t  2  p 

respectively,  with  r  the  position  vector  of  a  free  surface  fluid  particle,  g  the  acceleration  due  to 
gravity,  z  the  vertical  coordinate  (positive  upwards,  and  z  =  0  at  the  undisturbed  free  surface), 
pa  the  atmospheric  pressure,  p  the  fluid  density,  and  the  material  derivative  being  defined  as, 


D 

m 


d_ 

di 


+  u  ■  V 


4Note  that  such  a  redistribution  of  nodes  does  not  constitute  smoothing  (as,  e.g.,  used  in  LHC). 
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Waves  can  be  generated  by  simulating  a  piston  wavemaker  motion  on  the  “open  sea” 
boundary  of  the  computational  domain,  rri(i).  In  this  case,  wavemaker  motion  and  velocity 
[xp,  up]  are  specified  over  the  wavemaker  paddle  as, 

x  =  xp  ;  V</>  •  n  =  —  =  —Up  on  rVi(i)  (6) 

on 

where  overlines  denote  specified  values  (see  GS  for  detail). 

Along  the  stationary  bottom  T^,  and  on  other  fixed  boundaries  Yt2,  a  no-flow  condition  is 
prescribed  as, 

^  =  0  on  Tb  and  Vr2  (7) 


2.2  Boundary  discretization 

The  BIE  (2)  for  (</>,  |^)  is  solved  by  a  Boundary  Element  Method  (BEM,  Brebbia  (1978))  using 
a  set  of  Nr  collocation  nodes  on  the  boundary  and  Mr  =  M  +  Mf  higher-order  elements  to 
interpolate  in  between  the  collocation  nodes  (with  Mf,  the  number  of  elements  on  the  free 
surface,  and  M,  the  number  of  elements  on  the  rest  of  the  boundary).  In  these  boundary 
elements,  both  the  boundary  geometry  and  the  field  variables  (<j>,  |^)  are  interpolated  based  on 
nodal  values  and  on  polynomial  shape  functions  and/or  on  cubic  splines. 

Integrals  in  (2)  are  transformed  into  a  sum  of  integrals  over  the  elements  which  are  numer¬ 
ically  calculated  within  a  single  reference  element  T^,  with  intrinsic  coordinate,  £  £  [-1,+1],  A 
curvilinear  change  of  variable  is  defined  to  map  each  element  of  boundary  segment  Yk  onto  the 
reference  element.  We  hence  obtain  the  following  discretized  forms  for  the  integrals, 

/,  J^(v)G(x,xi)dV(x)  = 

Jr(x)  on 


N-p  M+Mf 


dsk  „  d</> 


E(E  Jrm)G(x(o,Xl)^-u)dc} 


j= 1  k= 1 

Nr  M+Mf  o  ,  Nr 

E(  E  4 r-t  =  ^K‘ 


d( 


dnyXj/ 


j=  1  k=l 


3=  1 


dlj'dn 


(8) 


dG(x,  xi) 
dn 


dT(x) 


m) 


dG(x((),  xi)  dsk 
]fn  ~di 


(0  d( 


Nr  M+M} 


Et  E 


Np 

E  K^i 


3=  1 


(9) 


in  which,  l  =  1, . . . ,  Nr’,  j  denotes  nodal  values  with  Nj(()  the  reference  element  shape  function 
at  node  j ;  and  dsk/d£  is  the  Jacobian  of  the  transformation  for  the  k- th  boundary  element. 
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Using  (8), (9),  the  discretized  form  of  the  BIE  (2)  reads, 


<>*  =  £  {K^-K^M  (10) 

0  =  i 

in  which,  l  =  l, ...  ,Nr.  Boundary  conditions  are  introduced  in  (10)  for  (f>  (see  time  stepping 
section)  and  for  (see  equations  (6),  (7)),  and  the  final  algebraic  system  is  assembled  by 
moving  unknowns  to  the  left  hand  side  and  keeping  specified  terms  in  the  right  hand  side.  The 
solution  of  this  system  is  simply  calculated  using  Kaletsky’s  LU  elimination  method.  Note  that 
coefficients  cti  are  computed  together  with  the  diagonal  coefficients  of  Knij  using  the  “rigid 
mode”  method  (Brebbia  (1978)).  This  method  has  been  shown  to  improve  both  the  conditioning 
and  the  accuracy  of  the  solution  of  the  resulting  algebraic  system  (GSS). 

Expressions  for  both  the  shape  functions  and  the  Jacobian  in  (8), (9)  are  given  in  GSS  for 
isoparametric  elements.  In  the  applications,  these  will  be  used  for  the  discretization  on  lateral 
boundaries  Tr| ,  Tr2,  and  on  the  bottom  Tf,  (Fig.  1).  Expressions  for  both  the  shape  functions  and 
the  Jacobian  on  the  free  surface  T  f  are  detailed  in  the  following  for  three  different  interpolation 
methods  (QS,  Mil,  and  MCI). 


(i)  Quasi- Spline  elements  (QS) 


As  discussed  in  the  introduction,  QS  elements  were  introduced  in  GS  to  improve  the  accuracy 
of  computations  in  highly  curved  regions  of  the  free  surface.  In  the  QS  elements,  shape 
functions  are  linear  polynomials  in  £  and  the  geometry  is  modeled  by  cubic  splines  (ensuring 
continuity  of  the  inter-element  slope).  To  represent  multiple-valued  free  surface  elevations 
occurring  in  overturning  waves,  the  geometry  is  defined  based  on  two  standard  single-valued 
spline  approximations,  for  [x(rj),  Tj]  and  [z(rj),  Tj ],  in  which  r  denotes  the  node  index, 


tOj)  =  U  =  and  £(t)  =  2(t  -  Tj)  -  1  (11) 


with  Nf  the  number  of  nodes  on  the  free  surface. 

For  the  k-\h  QS  element,  geometry,  xk(r)  =  ( xk(r ),  zk(r)),  is  thus  defined  as, 


xk(r)  =  (k  +  1  -  r)xk  +  (t  -  k)xk+l 


d2xk+ 1 
dr 2 


(12) 


A(t)  =  k  +  1  -  t  ;  B(t)  =  t  -  k;  t(£)  =  k  +  -  ^  (13) 

in  which  Tj  =  k  for  element  Tk  (k  =  1, . . . ,  MQ,  (xk,  xk+])  denote  element  k  nodal  coordinates, 
and  second-order  derivatives  of  these  are  obtained  from  the  spline  analyses.  End  node  slopes, 
dx / dr  and  dz/ dr,  required  as  boundary  conditions  in  the  spline  analyses,  are  in  general  specified 
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for  both  extremities  on  the  free  surface  based  on  4-node  cubic  polynomial  approximations.  In 
some  specific  cases,  however,  the  slope  at  one  extremity  of  the  free  surface  can  be  obtained 
from  the  physics  of  the  problem  (e.g.,  zero  slope  at  a  vertical  wall)  and  hence  can  be  explicitly 
specified  as  a  boundary  condition  in  the  spline  analyses. 


The  Jacobian  of  the  transformation  from  the  k- th  free  surface  element  to  the  reference 
element  is  defined,  using  (11)  as, 


dsk  dsk  9t 

lr(f)  =  ^(4(f) 

(14) 

(15) 

or, 

Qa k 

—  (£) 
d£ 

which  can  be  calculated  using  (12),  (13).  The  outward  unit  normal  vector  along  the  k- th  element 
is  similarly  obtained  as, 

„k 


lr  dzk  ,  dxk  ,  ,  dsk 


(16) 


(ii)  Cubic  Mid-Interval  Interpolation  method  (Mil) 

OSG  developed  a  full  cubic  interpolation  method,  the  Mil  method,  in  which  both  geometry  and 
field  variables  are  approximated  between  each  pair  of  nodes  (k,  k  +  1)  on  the  free  surface — 
except  for  the  first  and  last  pairs  including  corners — using  the  mid-section  of  a  cubic  (four-node) 
“sliding”  isoparametric  element.  At  both  extremities  on  the  free  surface,  the  first  two  nodes 
and  the  last  two  nodes  of  the  sliding  element  are  used  for  the  interpolation  in  between  the  first 
and  last  pairs  of  nodes,  respectively.  After  each  interpolation,  the  four-node  element  is  “slid” 
forward  along  the  free  surface,  which  ensures  local  continuity  of  the  inter-element  slope. 

Shape  functions,  to  be  used  in  (8),  (9),  are  cubic  polynomials  defined  in  the  sliding  element 
of  intrinsic  coordinate  fi  by, 

Nj(m)  =  Sij  with  m  =  (2 i  -  5)/3  (17) 

in  which  is  the  Kronecker’s  symbol;  i,j  =  1, . . . ,  4;  and  /il  are  intrinsic  nodal  coordinates 
for  the  sliding  element.  Hence, 

NM  =  4(1  _,,)(!)/ -1)  ;  JV2W= -4(1 -/)(!- 3^) 

NM  =  4(l_^)(1  +  3jl)  .  JV4(„)  =  4(1  + „)((>„*_!)  (18) 

The  mid- interval  segment  of  the  sliding  element  is  mapped  onto  the  usual  reference  element  T^ 

by, 

KO  =  M  +  (Hr  ~  (19) 
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in  which  (f ii,(ir )  denote  intrinsic  coordinates  for  the  leftward  and  rightward  nodes  of  the  Mil 
interval,  i.e.,  by  (17) :  (  — j,  +  j)  for  the  mid-interval,  and  (  —  1,  —  |)  or  (+|,  +1)  for  the  leftward 
and  rightward  extremities. 

The  Jacobian  of  the  transformation  from  the  k- th  free  surface  element  (i.e.,  segment  (k,  k  + 
1))  to  the  reference  element  is, 


dsk 

~dt 


(0 


,dzk 
'  dfi 


(KO)) 


2li 


(20) 


where,  using  (19),  d[i/d(  =  (fiT  —  hi)/ 2  =  1/3,  and, 


**(/*)  =  YlxiNM 

3= 1 


and 


(21) 


in  which  (xk,  j  =  1, . . . ,  4)  denote  coordinates  for  the  nodes  of  the  sliding  element  centered  on 
free  surface  segment  k,  and  N'-{n)  are  derivatives  of  shape  functions  (18). 


(iii)  Mixed  Cubic  Interpolation  method  (MCI) 

In  the  present  paper,  the  QS  and  the  Mil  methods  are  combined  into  the  “Mixed  Cubic  Interpo¬ 
lation  method”  (MCI)  for  which  field  variables  are  interpolated,  as  in  the  Mil  method,  using  the 
mid-section  of  a  cubic  “sliding”  isoparametric  element  and  the  geometry  is  modeled,  as  in  the  QS 
method,  using  two  cubic  spline  approximations,  for  [x{rj),  Tj\  and  [z{rj),  Tj\  (j  =  1, . . . ,  Nf). 
In  the  MCI  method,  equations  (1 1)  to  (16)  are  thus  used  for  the  modeling  of  free  surface  geom¬ 
etry  and  shape  functions  defined  by  (17)  to  (19),  are  used  for  interpolating  the  field 

variables. 

Both  the  accuracy  of  the  MCI  method  and  reasons  for  introducing  this  method  will  be 
discussed  in  the  applications. 


2.3  Numerical  integrations 

Non-singular  integrals  in  (8),(9)  (Ik  and  Ik  for  l  ^  j )  are  calculated  over  each  element 
k  using  a  standard  Gauss  quadrature,  usually  with  ten  integration  points  in  order  to  obtain 
sufficient  accuracy  for  higher-order  elements,  in  domains  with  large  aspect  ratios  (e.g.,  GS). 
Integrals  with  a  weak  logarithmic  singularity  (Ik  for  l  =  j )  are  first  transformed  into  a  standard 
form  by  isolating  the  singular  kernel  which  is  then  integrated  using  a  numerical  quadrature  exact 
for  the  logarithmic  singularity  (GSS).  Strong  singularities  in  the  integrals  (Ik  for  l  =  j )  are 

<’3 

removed  during  the  derivation  of  the  BIE  (2),  to  constitute  part  of  the  values  of  coefficients 
a(xi).  Although  non-singular,  resulting  integrals  may  still  have  a  highly  varying  kernel  when 
distance  |  x  —  xt  \  gets  small.  Accuracy  of  these  integrals  is  improved  by  first  performing  an 


9 


integration  by  parts  and  then  using  a  standard  Gauss  quadrature  rule  to  evaluate  the  remaining 
integral  (GSS). 

Quasi-singular  integrals  may  occur  near  comers  and  in  other  regions  of  the  free  surface,  like 
overturning  breaker  jets,  where  nodes  are  close  to  elements  on  different  parts  of  the  boundary. 
Although  non-singular,  these  integrals,  due  to  the  short  distance  from,  say,  collocation  node 
xi  to  the  boundary,  have  a  highly  varying  kernel  that  a  standard  Gauss  quadrature  rule  fails 
to  accurately  integrate  unless  a  prohibitive  number  of  integration  points  is  used  (GS).  Grilli 
and  Subramanya  (1994)  showed  that  the  loss  of  accuracy  of  Gauss  integrations  (with  ten 
integration  points)  for  such  quasi-singular  integrals  may  be  of  several  orders  of  magnitudes 
when  the  distance  to  the  collocation  node  becomes  very  small.  GS  developed  an  adaptive 
integration  scheme  based  on  a  binary  subdivision  of  the  reference  element  and  obtained  almost 
arbitrary  accuracy  for  the  quasi-singular  integrals,  when  increasing  the  number  of  subdivisions. 
This  method,  however,  can  be  computationally  expensive  and  Grilli  and  Subramanya  (1994) 
developed  a  more  efficient  method  that  essentially  redistributes  integration  points  around  the 
location  of  the  quasi-singularity  (point  of  minimum  distance  from  an  element  k  to  the  nearest 
collocation  node,  xi).  The  latter  method  will  be  used  in  the  present  applications  for  calculating 
quasi-singular  integrals. 


2.4  Node  regridding 

As  discussed  in  the  introduction,  free  surface  discretization  nodes  represent  fluid  particles  and, 
hence,  drift  in  the  direction  of  the  mean  mass  flux  of  the  flow,  thereby  affecting  resolution  of 
the  discretization.  To  either  add  and  redistribute  nodes  in  regions  of  poor  resolution  of  the 
free  surface  or  to  remove  and  redistribute  nodes  in  regions  of  node  concentration,  a  regridding 
technique  is  implemented  in  the  model  in  combination  with  the  MCI  interpolation  method.  This 
technique  redistributes  nodes  within  a  specified  boundary  section  based  on  a  constant  arc  length 
interval.  Field  variables  are  then  re-interpolated  at  the  new  locations  of  nodes  in  the  section. 

More  specifically,  a  node  arc  length  vector,  si  (with  l  =  1, ...,  Nf)  is  first  calculated  for  the 
old  locations  of  nodes,  xi ,  by  computing,  f( ds,  over  the  entire  free  surface  as  a  summation 
over  free  surface  elements.  A  boundary  section  is  then  defined  over  which  regridding  is  to  be 
performed  and  the  (arbitrary)  number  of  nodes  to  be  added  to  or  removed  from  this  section  is 
used  for  computing  the  new  arc  length  increment  within  the  section.  Based  on  this  increment, 
a  new  arc  length  vector,  sj,  is  calculated  forregridded  nodes  and  mapped  node  by  node  onto  s{ 
to  isolate  elements  k  (old  nodes  (k,  k  +  1))  within  which  new  nodes  are  located.  For  each  such 
element,  an  iterative  procedure  is  used  to  calculate  the  intrinsic  coordinate  £  corresponding  to 
a  given  component  of  s\  and,  based  on  this,  nodal  coordinates  xt  =  xk{()  are  re-calculated 
for  the  new  node  locations,  using  (12)  and  (13)  for  the  MCI  method,  and  similarly  for  the  field 
variables,  using  the  element  shape  functions  (18)  with  (19). 

Computations  are  pursued  as  before  in  the  regridded  discretization.  We  will  see,  however, 
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that  the  optimal  size  of  the  time  step  in  the  time  updating  part  of  the  model  is  function  of 
the  minimum  distance  between  free  surface  nodes.  Since  regridding  changes  this  distance, 
adjustments  must  also  be  made  to  the  time  step  in  order  for  the  computations  to  stay  accurate 
and  stable.  This  will  be  discussed  in  the  applications. 

As  mentioned  above,  discretization  nodes  gather  in  regions  of  flow  convergence.  When 
nodes  are  added  to  a  free  surface  section  and  regridded,  the  mean  distance  between  nodes 
in  this  section  is  decreased  and  quasi-singular  situations  are  likely  to  occur  even  more  often 
than  prior  to  regridding.  To  prevent  free  surface  nodes  from  getting  too  close  to  each  other 
(in  fact,  with  or  without  regridding)  and  to  limit  the  occurrence  of  such  quasi-singularities,  an 
adaptive  regridding  procedure  is  implemented  in  the  model.  When  the  distance  between  two 
nodes  becomes  50%  smaller  (or  larger)  than  the  distance  between  neighboring  nodes,  nodes  are 
locally  regridded  to  make  these  distances  identical.  This  procedure  will  also  be  discussed  and 
used  in  the  applications. 


2.5  Time  stepping  and  global  accuracy 


Free  surface  boundary  conditions  (3)  and  (4)  are  integrated  at  time  t  to  establish  both  the  new 
position  and  the  boundary  conditions  on  the  free  surface  at  subsequent  time  (t  +  At)  (with 
At  the  time  step).  Second-order  Taylor  series  expansions  are  used  to  express  both  the  new 
position  r(t  +  At)  and  the  potential  4>(r(t  +  At))  on  the  free  surface,  in  an  Eulerian-Lagrangian 
formulation  based  on  the  material  derivative  (5), 

Dr  (At)2  D2r 

r(i  +  At)  =  r(i)  +  At—  (<)  +  +  0[(Ai)3]  (22) 

?(r(t  +  At),  t  +  At)  =  m  +  Aij^(i)  +  AAtAtf)  +  0[(  A  i)3]  (23) 

Coefficients  in  the  Taylor  series  are  expressed  as  functions  of  the  potential,  its  time  derivative,  and 
the  normal  and  tangential  derivatives  of  both  of  these  along  the  free  surface  (GSS;  Grilli  (1993)). 
This  method  has  the  advantages,  compared  to  Runge-Kutta  (RK)  and  Adams-Bashforth-Moulton 
(ABM)  time  stepping  schemes  used  in  earlier  models  (LHC;  VB;  NMP),  of  both  being  explicit 
and  accounting  for  spatial  derivatives  along  the  free  surface  (d/ds,  d2/ds2)  in  the  calculation 
of  values  at  ( t  +  At).  This  provides  a  better  stability  of  the  computed  solution  and  makes  it 
possible  to  use  larger  time  steps,  for  a  similar  accuracy,  thus  also  making  the  overall  solution 
more  efficient  (see  also  discussions  in  RO,  OSG). 


To  calculate  coefficients  in  series  (22),  (23),  normal  and  temporal  derivatives  of  the  potential 


are  obtained  from  the  solution  of  two  BIE’s  of  the  form  (2),  for  (</>,  |^)  and  (?£,  -2-2-),  and 


'  dcj>  dl(f> 


tangential  derivatives  are  computed  using  a  five-node  4th-order  “sliding”  element,  independent 
from  the  BEM  discretization.  Boundary  conditions  for  the  second  BIE  for  are  obtained 


from  the  solution  of  the  first  BIE  for  (f>.  Note  that  both  BIE’s  correspond  to  the  same  boundary 


geometry  and,  hence,  have  the  same  discretized  form  (10).  Therefore,  the  solution  of  the  second 
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BIE  only  takes  a  small  fraction  of  the  time  (typically  a  few  percent)  needed  for  the  solution  of 
the  first  BIE  which,  again,  makes  the  present  time  stepping  method  very  efficient,  particularly, 
compared  to  higher-order  RK  or  ABM  schemes  which  often  require  multiple  evaluations  of  the 
BIE  (2)  for  several  intermediate  times  per  time  step. 

In  the  applications,  global  accuracy  of  the  computations  is  checked  at  each  time  step  by 
computing  errors  in  total  volume  m  and  energy  e  for  the  generated  wave  train.  Numerical  errors 
are  function  of  the  size  (i.e.,  distance  between  nodes)  and  degree  (i.e.,  quadratic,  cubic,...)  of  the 
boundary  elements  used  in  the  spatial  discretization,  both  of  which  control  the  accuracy  of  the 
BEM  solution  of  Laplace’s  equation,  and  of  the  size  of  the  selected  time  step,  which  controls  the 
accuracy  of  the  time  stepping.  GS  developed  a  criterion  for  selecting  the  optimum  time  step  size 
in  the  model  based  on  results  of  computations  made  in  various  spatio-temporal  discretizations  for 
the  propagation  over  a  given  time  (i.e.,  distance)  of  a  large  solitary  wave.  Using  QS  elements  on 
the  free  surface  and  quadratic  isoparametric  elements  elsewhere  on  the  boundary,  they  showed 
that,  for  a  constant  time  step  A ta,  errors  in  m  and  e  reach  a  minimum  when  the  mesh  Courant 
number,  CQ  =  \fgh  Ata/ Axa  ~  0.5  (in  which  Ax0  denotes  the  initial  distance  between  nodes 
on  the  free  surface  and  h  the  characteristic  depth).  Thus,  despite  more  accurate  results  in  the 
time  updating  (22), (23)  when  a  smaller  time  step  is  used  (truncation  errors  are  0[(At 0)3]),  the 
model  may  provide  larger  overall  errors,  in  a  given  initial  discretization,  when  the  time  step  is 
too  small.  The  reason  for  this  is  that,  for  very  small  time  steps,  too  many  computational  steps 
are  needed  for  the  wave  to  cover  the  specified  time  (or  distance)  of  propagation  in  the  model, 
which  leads  to  an  “accumulation”  of  numerical  errors  (i.e.,  truncation,  round-off)  resulting  from 
the  solution  of  Laplace’s  equation.  Hence,  in  such  cases,  overall  accuracy  can  only  be  improved 
by  reducing  Axa,  thus  making  the  solution  of  Laplace’s  equation  more  accurate,  and  achieving 
optimum  Ca . 


Based  on  these  results,  GS  proposed  an  adaptive  time  step  method  applicable  to  highly 
transient  waves — like  breakers — for  which  the  distance  between  nodes  widely  varies.  In  this 
method,  time  step  At  is  adaptively  calculated  as  a  function  of  time  t  based  on  the  optimum 
mesh  Courant  number,  Ca,  and  on  the  instantaneous  minimum  distance  between  nodes  on  the 
free  surface,  A  |  r{t)  |mm, 


At 


min 

•Jgh 


(24) 


Similar  calculations  will  be  presented  in  the  first  application  in  this  paper,  using  the  Mil  or 
the  MCI  method  for  the  interpolation  on  the  free  surface,  in  order  to  determine  the  corresponding 
optimum  value  of  Ca,  if  any. 


2.6  Discretized  boundary  conditions  at  corners 

In  the  following,  we  assume  that  initial  boundary  conditions  are  selected  in  the  model  in  such 
a  way  that  the  flow  at  comers  does  not  have  an  initial  (mathematical)  singularity  (see  Section 
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3.3  for  details).  As  discussed  in  the  introduction,  well-posed  discretized  boundary  conditions 
must  be  specified  at  corners  to  avoid  occurrence  of  numerical  singularities  during  the  course  of 
computations  5. 

Boundary  conditions  and  normal  directions  are  in  general  different  for  each  side  of  comers 
in  the  computational  domain  (e.g.,  on  a  wavemaker  and  on  the  free  surface;  see  Fig.  1).  To  be 
able  to  specify  such  differences,  corners  are  represented  by  double-nodes  for  which  coordinates 
are  identical  but  normal  vectors  are  different  (GSS;  GS).  Thus,  two  different  discretized  BIE’s 
(eq.  (10))  are  expressed  for  each  node  of  a  corner  double-node,  say  :  (i)  for  l  =  p,  on  the 
wavemaker  paddle;  and  (ii)  for  l  =  /,  on  the  free  surface.  Since  potential  must  be  unique  at  a 
given  location,  however,  one  of  these  two  BIE’s  must  be  modified  in  the  final  discretized  system 
to  explicitly  satisfy,  4>p  =  4>f  (i.e.,  “continuity  of  the  potential”,  GS). 

In  addition,  for  free  surface  corners,  the  velocity  vector  must  also  be  unique,  i.e.  up  =  Uf, 
for  both  corner  (double)  nodes  to  move  to  identical  locations  through  time  updating.  GS  and 
OSG  showed  that,  if  this  “velocity  continuity”  is  not  explicitly  specified  in  the  discretized  BIE 
system  (10),  large  numerical  errors  will  occur  at  corners  which  will  grow  even  larger  through 
time  stepping  and,  eventually,  will  lead  to  instability  of  the  solution,  particularly  close  to  a 
strongly  moving  lateral  boundary  6.  Hence,  they  derived  discretized  relationships  expressing 
velocity  uniqueness  (compatibility)  at  corners  for  all  cases  of  mixed  Dirichlet-Neuman  boundary 
conditions.  These  relationships,  in  fact,  make  the  representation  of  the  solution  compatible  (i.e., 
consistent)  on  both  sides  of  a  corner  in  the  sense  of  the  study  by  Gray  and  Lutz  (1990)  and 
effectively  eliminate  (numerical)  singularities  in  the  discretized  solution. 

When  using  compatibility  conditions  at  comers  for  solving  mixed  boundary  value  problems 
in  simple  rectangular  domains,  GS  showed  that  numerical  errors  at  corners  could  be  reduced  to 
almost  arbitrarily  small  values.  At  the  intersection  between  a  piston  wavemaker  and  the  free 
surface,  for  instance,  they  used  a  velocity  compatibility  condition  to  specify  a  (corrected)  value 
of  the  tangential  velocity,  as  a  function  of  both  the  normal  velocity,  obtained  from  the 

solution  of  the  BIE  (10)  at  each  time  step,  and  the  (known)  wavemaker  velocity,  ^  =  — up(t ) 
(see  (6)),  before  calculating  time  updating  of  the  free  surface.  This  condition  reads, 


ds 


dn 


tan  /3f 


dcj)p 

dn 


sec  (3  f 


(25) 


in  which  (3f  is  the  angle  between  the  free  surface  and  the  horizontal.  The  velocity  obtained  in 
(25)  is  thus  used  in  replacement  of  the  tangential  velocity  otherwise  calculated  in  the  4th-order 
sliding  polynomial,  to  make  velocity  components  a  posteriori  compatible  for  both  sides  of  the 
comer. 

5  Gray  and  Lutz  (1990)  and  Gray  and  Manne  (1993)  showed  that,  provided  the  discretization  of  field  variables  on 
both  sides  of  a  corner  is  consistent,  terms  which  are  potentially  singular  at  corners  will  cancel.  In  turn,  ill-posedness 
of  corner  boundary  conditions  through  inconsistent  discretization  will  create  the  equivalent  of  hypersingularities. 

6This  is  also  consistent  with  the  conclusions  in  Gray  and  Lutz  (1990)  and  Gray  and  Manne  (1993). 
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By  performing  additional  numerical  tests,  OSG  showed  that  using  (25)  a  posteriori  may 
not  be  accurate  enough  for  strongly  moving  lateral  boundaries.  Hence,  they  further  developed 
this  equation  to  a  priori  specify  velocity  compatibility  at  each  time  step  within  the  discretized 
BIE  system  (10),  before  the  solution  is  calculated.  To  do  so,  they  first  expressed  the  tangential 
velocity  in  (25)  as  a  linear  combination  of  nodal  values  of  the  potential  within  a  4th-order  sliding 
element  located  over  the  first  five  nodes  of  the  free  surface  (k  =1), 

Bch  5  3<*k  5 

-P  =  {T,*riMM/-grM  =  (26) 

j=  i  *=  j=  i 

in  which,  by  (17),  p\  =  —  1  for  the  first  node  in  the  element,  and  N'-{p)  denote  first  derivatives 
of  the  sliding  element  quartic  shape  functions,  obtained  by  (17)  with  i,  j  =  5.  Combining  (25) 
and  (26),  they  finally  expressed  the  corner  compatibility  equation  as, 

T\$P  ~  ^r-  tan/3/  =  up  sec/3/  -  (27) 

dn  j=- 

in  which  </>i  =  </>/  =  cf>p,  and  the  last  term  in  the  right  hand  side  is  a  linear  combination  of 
values  of  the  free  surface  potential  specified  using  time  updating  equation  (23).  Equation  (27) 
is  introduced  in  line  l  =  p  of  the  discretized  system  (10),  in  replacement  of  the  BIE  for  the 
wave  maker  side  of  the  comer,  and  the  other  comer  equation  in  line  l  =  f  of  the  system  is,  as 
before,  the  discretized  BIE  for  the  free  surface  side  of  the  comer.  Although  prescribed  on  the  free 
surface  by  the  time  updating  (23),  the  potential  in  the  corner,  </>p  =  </>/,  is  now  considered  as  an 
unknown  to  be  calculated  as  part  of  the  system  solution  in  order  to  satisfy  velocity  compatibility. 
Hence,  coefficients  of  4>p  and  </>/  are  lumped  together  for  lines  l  ^  p  of  the  system  (10). 

Compatibility  equations  similar  to  (27)  are  also  expressed  for  ^  and  for  all  the  other  cases 
of  mixed  boundary  conditions  at  corners  (Dirichlet  and  Neuman  types)  discussed  by  GS.  OSG 
showed  that  using  equations  such  as  (27)  at  corners  significantly  improves  both  the  accuracy  and 
the  stability  of  the  numerical  solution,  particularly  for  cases  of  waves  generated  by  wavemakers. 
These  new  extended  compatibility  conditions  have  been  implemented  in  the  present  model  and 
their  accuracy  is  compared  to  the  old  comer  compatibility  conditions  (type  (25)),  in  the  last 
application  in  the  next  section. 


3  Applications 

3.1  Global  accuracy  for  solitary  wave  propagation  over  constant  depth 

Computations  are  carried  out  as  in  GS  to  determine  the  optimal  mesh  Courant  number  Ca 
corresponding  to  minimum  numerical  errors  for  given  discretizations  and  interpolation  methods 
on  the  free  surface.  Numerically  exact  solitary  waves  obtained  using  the  fully  nonlinear  method 
by  Tanaka  (1986)  are  used  as  incident  waves  (see  GS  for  detail)  and  are  propagated  over 
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constant  depth  h  in  various  spatio-temporal  discretizations  (Fig.  2).  Such  solitary  waves 
should  keep  permanent  form,  constant  volume  and  total  energy  while  propagating  in  the  model. 
Computational  errors  with  respect  to  Tanaka’s  results  for  each  wave  thus  give  a  measure  of 
discretization  and  time  step  effects  on  global  numerical  accuracy,  for  different  interpolation 
methods  used  in  the  free  surface. 

Solitary  waves  of  height  H/h  =  0.3  to  0.6  are  used  in  the  computations  and  computational 
data  are  identical  to  those  for  the  same  applications  in  GS.  The  domain  aspect  ratio  is  L/h  =  28 
(with  h  =  1).  Three  different  spatial  discretizations  are  used  in  the  computations,  with  initial 
distance  between  nodes  on  the  free  surface,  AxJ,  =  Axa/h  =  0.15,  0.20,  0.25.  In  addition, 
A z'  =  0.25  on  the  lateral  boundaries  in  each  discretization  and  Ax'  =  0.25  on  the  bottom,  in 
the  first  set  of  results,  and  0.40  in  the  next  two  ones.  The  QS,  Mil,  and  MCI  cubic  interpolation 
methods  are  successively  used  on  the  free  surface,  while  three-node  quadratic  isoparametric 
elements  are  used  on  other  boundaries  of  the  domain  (Fig.  2).  Ten  Gauss  points  are  used  per 
element  and  adaptive  integration  is  specified  in  comer  elements. 

Figs.  3a-c  show  results  for  a  solitary  wave  of  height  H/h  =  0.50.  The  incident  wave  is 
introduced  in  the  model  with  its  crest  at  x'  =  14  and  is  then  propagated  over  5  dimensionless 
time  units  (Fig.  2).  Maximum  absolute  errors  in  wave  volume  m,  total  energy  e,  and  wave  height 
H  are  calculated  over  this  propagation  time  with  respect  to  Tanaka’s  solution  (m  =  1.7914787, 
e  =  0.6157121,  and  H  =  0.5000000).  Four  different  constant  time  steps  are  used  in  the 
computations7,  A t'0  =  Atayjg/h  =  0.025,  0.05,  0.10,  and  0.20.  Twelve  cases  with  different 
spatio-temporal  discretization  are  thus  calculated  using  each  interpolation  method  and  errors 
are  plotted  in  the  figures,  in  negative  logarithmic  scales,  as  a  function  of  Ax'a  (symbols),  the 
interpolation  method  (lines),  and  Ca  =  At'0/ Ax'a.  Note  that  results  with  the  QS  method  are 
identical  to  those  in  GS. 

Results  in  Figs.  3a-c  show  that  the  Mil  and  the  MCI  methods  give  similar  errors  in  most 
cases  (except  for  the  largest  Courant  number)  and  that,  for  CQ  <  0.5,  errors  are  reduced  by  up  to 
3  orders  of  magnitude  as  compared  to  the  QS  method.  It  should  be  stressed  that  these  reductions 
in  errors  are  obtained  for  essentially  the  same  computation  time.  With  the  QS  method,  as 
expected  from  both  the  conclusions  in  GS  and  the  discussion  in  Section  2.5,  errors  are  minimum 
or  level  up  for  Ca  ~  0.50.  With  the  Mil  and  the  MCI  methods,  only  errors  in  volume  in  the 
coarser  discretization  reach  a  minimum  (Fig.  3a)  and,  for  the  other  tested  cases,  errors  either 
do  not  reach  a  minimum  or  more  or  less  level  up  for  Ca  ~  0.20-0.35.  This  is  because  the 
more  accurate  solution  of  Laplace’s  equation  with  the  MCI  and  Mil  methods  leads  to  a  slower 
rate  of  error  accumulation,  when  many  time  steps  are  used.  To  avoid  using  too  many  small 
computational  time  steps,  however,  the  upper  bound  of  this  interval,  i.e.,  Ca  ~  0.35,  will  be 
used  in  later  applications  of  the  Mil  and  MCI  methods,  in  combination  with  the  adaptive  time 
step  procedure  (24).  As  pointed  out  in  GS  and  seen  in  Fig.  3,  accuracy  can  only  be  significantly 

7Note  that,  as  discussed  in  Section  2.5,  the  smaller  the  time  step,  the  larger  the  number  of  time  steps  necessary 
to  propagate  the  wave  over  the  given  time. 
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improved  around  this  selected  value  of  CD  by  using  a  smaller  initial  Ax'a. 

The  influence  of  initial  wave  height — a  measure  of  nonlinearity — on  model  accuracy  is 
tested  next  and  results  are  plotted  in  Fig.  3d.  Solitary  waves  of  height,  H/h  =  0. 3-0.6,  are 
propagated  in  the  same  domain  as  above  (Fig.  2),  using  Ax'0  =  0.15,  CQ  =  0.35,  the  adaptive 
time  step  procedure  (24),  and  the  MCI  method  for  the  interpolation  on  the  free  surface.  Fig.  3d 
gives  the  same  types  of  errors  as  in  Figs.  3a-c.  One  can  see  that  the  model  is  almost  uniformly 
accurate  over  the  whole  range  of  tested  wave  heights,  up  to  fairly  steep  waves.  One  can  also 
see,  comparing  results  for  H/h  =  0.5,  that  errors  are  slightly  smaller  when  using  the  adaptive 
time  step  scheme  (Fig.  3d)  than  when  using  fixed  time  steps  (Figs.  3a-c). 


3.2  Solitary  wave  shoaling  and  breaking  over  a  gentle  slope 

Solitary  waves  are  often  used  as  a  simple  model  for  studying  propagation,  shoaling,  and  breaking 
of  extreme  waves  in  shallow  water.  Solitary  waves  closely  model  tsunamis  and  can  also  be 
used  to  represent  surf-zone  waves  (see  GSSV  and  Grilli,  Svendsen  and  Subramanya  (1994)  for 
detail). 

Computations  of  shoaling  and  runup  of  fully  nonlinear  solitary  waves  over  a  slope  were 
reported  by  Svendsen  and  Grilli  (1990)  using  the  QS  method  and  by  OSG,  using  the  Mil 
method.  These  computations,  however,  were  not  performed  for  very  gentle  slopes  which  are 
more  demanding  cases,  due  to  the  quasi-singular  situations  in  the  BEM  resulting  from  the  very 
narrow  geometry  in  the  upper  part  of  the  slope.  Grilli  and  Subramanya  ( 1 994)  developed  efficient 
methods  for  calculating  quasi-singular  integrals  and,  using  the  MCI  method  for  the  interpolation 
on  the  free  surface,  reported  quite  accurate  results  for  the  shoaling  of  solitary  waves  over  a  1:35 
slope.  GSSV,  following  the  same  approach,  calculated  many  detailed  features  of  the  flow  under 
the  crest  of  shoaling  solitary  waves,  up  to  initiation  of  breaking  by  overturning  of  the  crest. 
In  all  these  cases,  termination  of  calculations  always  occurred  when,  due  to  flow  convergence 
in  poorly  resolved  breaker  jets,  two  nodes  moved  very  close  to  each  other  thus  creating  an 
almost  singular  situation  in  the  BEM,  which  led  to  very  large  errors  and  to  the  instability  of 
computations.  In  fact,  in  all  these  applications,  computations  had  to  be  stopped  a  short  time  after 
the  wave  crest  started  overturning,  and  no  detailed  analysis  of  the  breaker  jet  could  be  made  8. 

Errors  in  the  BEM  solution  can  be  reduced  in  breaker  jets  without  using  a  finer  initial 
discretization  over  the  whole  free  surface,  which  would  be  too  computationally  expensive,  by 
:  (i)  using  the  more  accurate  MCI  elements  in  the  discretization  instead  of  QS  elements;  (ii) 
improving  the  resolution  by  selectively  adding  nodes  to  the  wave  crest  region  (i.e.,  to  only  a 
small  portion  of  the  free  surface);  and  (iii)  limiting  the  occurrence  of  quasi-singular  situations 

8  In  these  computations,  the  poor  resolution  in  the  breaker  jets  resulted  from  the  large  horizontal  extension  of 
the  discretized  free  surface  in  the  computational  domain  (sometimes  more  than  100  water  depths)  which  prevented 
a  sufficient  number  of  nodes  from  gathering  in  the  wave  crest  region,  before  overturning  was  initiated,  unless  a 
prohibitive  number  of  nodes  was  initially  used  over  the  whole  free  surface. 
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by  preventing  nodes  from  moving  too  close  to  each  other  (adaptive  regridding).  Benefits  of 
using  methods  (i)  to  (iii)  are  illustrated  in  the  following  computations  for  the  shoaling  of  solitary 
waves  of  initial  height,  H'a  =  H0/h0  =  0.20  and  0.40,  over  a  1:35  slope  (i.e.,  in  a  set-up 
similar  to  the  applications  in  GSSV).  For  each  wave,  two  cases  are  calculated,  either  using  the 
QS  method  or  the  MCI  method  for  the  interpolation  on  the  free  surface,  and  their  accuracy  is 
compared.  Node  regridding  techniques  introduced  in  Section  2.4  are  then  used  to  increase  the 
resolution  of  the  discretization  in  breaker  jets,  while  preventing  nodes  from  moving  too  close  to 
each  other.  This  allows  calculations  to  be  carried  out  for  a  longer  time,  within  the  same  overall 
accuracy  9  (e.g.,  1%  maximum  error  on  wave  mass). 


(i)  Comparison  between  the  QS  and  MCI  methods 

The  computational  domain  is  sketched  in  Fig.  1.  Since  all  tested  waves  break  before  reaching 
the  top  of  the  slope,  the  upper  part  of  the  slope,  where  the  domain  geometry  becomes  very 
narrow,  was  replaced  by  a  small  shelf  of  depth  h  =  0.1  ha  (i.e.,  for  x'  >  41.5).  Unnecessary 
refinements  of  the  discretization  and  of  the  integration  methods  in  the  upper  part  of  the  slope 
are  thus  avoided.  It  can  be  shown  that  reported  results  are  not  influenced  by  this  slight  change 
in  geometry. 

There  are  Nr  =  429  nodes  in  the  discretization,  including  Nf  =  226  nodes  on  the  free 
surface  (with  initial  spacing  Ax'a  =  0.20)  representing  Mf  =  225  QS  or  MCI  elements.  There 
are  M  =  100  quadratic  isoparametric  elements  on  the  bottom  and  lateral  boundaries.  The 
distance  between  nodes  on  the  bottom  is  0.5  in  the  constant  depth  region  ( x '  <  10)  and  this 
distance  is  progressively  reduced  over  the  slope  to  0.40,  0.25,  0.20,  0.15  and  0.10,  in  order  to 
get  increased  resolution  when  depth  decreases.  On  the  shelf  bottom,  the  distance  between  nodes 
is  0.15.  Quasi-singular  integration  methods  by  Grilli  and  Subramanya  (1994)  are  used  in  all 
comer  elements  and  in  free  surface  and  bottom  elements  for  which  the  local  depth  is  equal  to  or 
smaller  than  the  element  length  on  the  free  surface,  Ax'a  (i.e.,  for  x'  >  38).  The  mesh  Courant 
number  is  set  to  CQ  =  0.35  for  both  interpolation  methods,  which  corresponds  to  A t'Q  =  0.070. 
With  these  data,  the  CPU  time  per  time  step  is  12.9  sec  (IBM9000/3)  for  both  the  QS  and  the 
MCI  methods  10 . 

Incident  solitary  waves  are  generated  :  (i)  using  the  numerical  piston  wavemaker  on  the 
leftward  lateral  boundary  of  the  domain,  for  H'a  =  0.20;  or  (ii)  using  the  exact  method  by 
Tanaka  (1986)  in  the  constant  depth  region  prior  to  the  slope,  for  H'a  =  0.40  (with  the  initial 
crest  at  x'  =  9).  Figs.  4a  (curves  a-d)  and  4c  (curves  A-D)  show  stages  of  wave  shoaling 
and  breaking  calculated  with  the  MCI  method  in  the  initial  discretization,  up  to  slightly  prior 

9  Grilli,  Svendsen  and  Subramanya  (1994)  and  Wei,  Kirby,  Grilli  and  Subramanya  (1995)  also  reported  similar 
calculations  of  flow  characteristics  of  shoaling  solitary  waves,  up  to  and  beyond  the  breaking  point,  over  slopes 
1:8  down  to  1:100.  These  references  can  be  consulted  for  more  results  and  discussions  on  the  physics  of  wave 
processes. 

10It  is  again  seen  that  the  more  accurate  MCI  method  is  not  more  computationally  expensive  than  the  QS  method. 
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the  instant  of  instability  of  the  computations.  The  breaking  point,  defined  as  the  crest  location 
for  which  the  free  surface  slope  becomes  vertical  in  the  wave  front,  is  located  at  x'b  =  36.3 
(reached  at  t'b  =  44.53)  and  x'b  =  30.4  (reached  at  tb  =  18.33),  for  each  wave  respectively. 
Corresponding  breaking  wave  heights  and  breaking  indices  are  :  (i )Hb  =  0.364,  Hb/hb  =  1.46; 
and  (ii)  Hb  =  0.587,  Hb/hb  =  1.41,  respectively.  During  propagation,  time  step  reduces  to 
At'  =  0.019  and  0.0003,  at  the  time  of  breaking,  for  each  wave  respectively.  The  total  number 
of  time  steps  up  to  the  breaking  point  is  950  and  750  and  the  average  time  step  is  0.047  and 
0.024  for  each  case,  respectively. 

Figs.  5a  and  5b  give  numerical  errors  in  wave  mass  and  total  energy  for  the  QS  and  the 
MCI  methods.  For  most  of  the  wave  propagation,  these  errors  are  much  smaller  using  the 
MCI  method  than  using  the  QS  method.  Hence,  for  a  pre-specified  maximum  numerical  error 
(e.g.,  1%),  wave  propagation  can  be  calculated  for  a  longer  time  with  the  MCI  method  before 
computations  are  deemed  inaccurate.  Fig.  5b  shows  that,  for  both  interpolation  methods,  the 
energy  error  first  reaches  the  1%  threshold  and,  as  expected,  this  occurs  quite  earlier  with  the 
less  accurate  QS  method.  Arrows  in  Fig.  5a  show  locations  for  points  of  maximum  energy 
error.  Fig.  5c  gives  relative  wave  heights  H/h  as  a  function  of  x' .  for  both  interpolation 
methods.  Arrows  again  indicate  locations  at  which  energy  errors  reach  1%  in  the  computations. 
As  expected,  using  the  MCI  method,  shoaling  is  calculated  for  a  longer  time  within  the  required 
accuracy11. 

In  the  last  stages  of  shoaling,  discretization  nodes  converge  towards  wave  crest  regions 
while  the  distance  between  nodes  increases  in  the  front  and  in  the  back  of  the  waves.  This  is  quite 
apparent  in  Figs.  4b  and  4d  (dashed  curves  d  and  D).  Despite  this  node  concentration,  however, 
when  breaker  jets  start  forming  at  the  tip  of  the  crests,  the  discretization  is  not  sufficiently 
resolved  to  allow  formation  of  large  scale  plunging  jets.  Hence,  numerical  errors  quickly 
increase  leading  to  instability  of  computations  (as  seen  in  Figs.  5a-b).  This  can  be  improved 
using  the  node  regridding  techniques. 


(ii)  Regridding  of  the  breaker  jet 

Node  regridding  is  used  to  improve  resolution  of  computations  (and  thus  numerical  accuracy) 
in  the  last  stages  of  wave  shoaling  and  breaking.  For  both  waves,  40  nodes  are  added  to  a 
free  surface  segment  containing  the  wave  crest,  at  a  time  slightly  before  reaching  the  breaking 
point  and  the  steep  increase  in  computational  errors  observed  in  Figs.  5a-b  (i.e.,  at  t'  =  44.12 
and  18.15  respectively).  These  nodes  are  regridded  to  constant  arc  length  intervals  together 
with  “old”  nodes  in  the  sections.  Regridded  sections  thus  span  from  x'  =  33.26  to  37.80 
and  x'  =  25.84  to  32.00,  and  the  number  of  “old”  nodes  in  these  sections  is  31  and  37,  for 

1 1  Note  that  (more  accurate)  shoaling  curves  calculated  with  the  MCI  method  are  slightly  lower  than  with  the  QS 
method.  This  agrees  to  within  2%  with  measurements  by  GSSV,  up  to  the  breaking  point. 
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each  wave  respectively  12 .  After  regridding,  the  total  number  of  nodes  in  the  discretization  is 
469  and  the  CPU  time  per  time  step  increases  to  14.8  sec.  To  reduce  quasi-singularities  due 
to  node  proximity  in  regridded  breaker  jets,  the  adaptive  regridding  procedure  is  used  for  the 
computations  in  the  regridded  discretizations  to  check  distances  between  nodes  and  regrid  nodes 
two-by-two  if  needed. 

Due  to  regridding,  the  minimum  distance  between  nodes  increases  in  the  discretizations. 
Hence,  assuming  that,  for  stability  reasons,  computations  are  resumed  in  the  regridded  domain 
with  a  time  step  equal  to  its  value  prior  to  regridding  (i.e.,  At'  =  0.0103  and  0.0014  for 
each  wave  respectively),  equation  (24)  implies  that  the  mesh  Courant  number  must  decrease. 
Although  computations  using  a  smaller  Courant  number  are  still  accurate  (see  Section  3.1),  this 
may  lead  to  an  unnecessary  large  number  of  time  steps  and,  eventually,  to  an  accumulation  of 
numerical  errors.  To  avoid  this  problem,  after  regridding,  the  Courant  number  is  increased  back 
to  0.35  over  10  computational  steps,  following  a  smooth  “tanh”  law  variation,  and  the  time 
step  is  correspondingly  increased.  After  regridding,  280  and  200  additional  computational  steps 
are  thus  calculated  for  each  wave,  respectively,  and  corresponding  time  step  values  reduce  to 
At'  =  0.0043  and  0.0061.  Wave  profiles  computed  in  the  regridded  discretizations  are  shown 
in  Fig.  4,  curves  e-k  and  E-I,  and  corresponding  errors  in  wave  mass  are  given  in  Fig.  5d  and 
compared  to  errors  obtained  for  the  same  waves  without  regridding  (note  that  errors  in  wave 
energy  are  about  5  times  as  large  as  these).  Last  plotted  profiles  in  Fig.  4  (curves  k  and  I) 
correspond  to  a  1%  error  in  wave  mass  13 . 

Fig.  5d  shows  that,  within  a  1%  accuracy  on  wave  mass,  computations  can  be  pursued 
for  a  longer  time  and  over  a  larger  distance  of  propagation  in  the  regridded  discretizations,  up 
to  the  theoretical  limit  of  FNPF  theory  (i.e.,  the  instant  the  tip  of  the  breaker  jet  is  about  to 
impact  the  free  surface  :  curves  k  and  I),  for  only  a  very  marginal  increase  in  computational 
effort.  Results  in  Figs.  4b  and  4d,  in  the  regridded  discretizations,  show  the  development  of 
well  resolved  breaker  jets.  Fig.  5c  shows  that,  in  this  post-breaking  region,  breaking  indices 
reach  an  almost  constant  value.  Both  waves,  finally,  exhibit  well  known  features  of  plunging 
breakers  over  gentle  slopes  (Peregrine  (1983)),  with  a  large  scale  plunging  of  the  wave  crest 
and  quasi-elliptical  areas  enclosed  below  the  breaking  jets  (“pipelines”).  Grilli,  Svendsen  and 
Subramanya  (1994)  discussed  the  physics  of  these  and  other  similar  results  and  showed  that 
they  agree  quite  well  with  experimental  data  from  various  sources. 

12Free  surface  shapes  and  node  locations  before  and  after  regridding  are  very  similar  to  those  shown  in  curves 
d,e  and  D.E  of  Figs.  4a-b,  for  each  wave  respectively. 

13Note  that  computations  by  Subramanya  and  Grilli  (1994),  in  which  regridding  was  performed  at  the  breaking 
point  using  a  larger  number  of  “new”  nodes  (49)  and  over  a  narrower  region  of  the  free  surface  (mostly  in  front 
of  the  wave  crest),  showed  that  both  numerical  errors  could  be  cut  to  about  one-fourth  the  present  case  and  free 
surface  profiles  obtained  were  closely  identical  to  those  in  Fig.  4.  This  illustrates  the  convergence  of  numerical 
procedures  used  for  computing  the  solution  in  breaker  jets. 
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3.3  Wave  breaking  induced  by  a  moving  boundary 

Problems  with  strong  interactions  between  a  moving  boundary  and  a  free  surface  are  quite 
frequent  and  important  ones  in  ocean  engineering.  To  name  but  a  few  of  these  :  the  generation 
of  a  (possibly  breaking)  bow  wave  by  a  forward  moving  vessel,  the  slamming  of  a  ship  hull  on 
the  ocean  surface,  or  the  generation  of  waves  by  a  wavemaker  in  a  laboratory  wave  tank. 

When  a  lateral  vertical  boundary  of  the  computational  domain,  say  Tr  |  in  Fig.  1,  is 
moved  from  a  state  of  rest  with  sufficient  acceleration,  the  free  surface  close  to  the  intersection 
with  the  boundary  quickly  rises  upward  and  also  acquires  a  large  horizontal  acceleration. 
Possible  mathematical  or  numerical  singularities  of  the  solution  in  the  comer  representing  the 
intersection  between  the  free  surface  and  the  lateral  boundary  were  discussed  in  the  introduction 
and  in  Section  2.6.  It  was  seen  that,  provided  initial  boundary  conditions  are  well-posed  (thus 
eliminating  initial  (mathematical)  singularities)  and  BEM  numerical  integrations  are  accurate, 
compatibility  conditions  could  be  used  to  ensure  instantaneous  well-posedness  of  boundary 
conditions  at  corners  and  avoid  the  occurrence  of  (numerical)  singularities  through  time  stepping. 
Doing  so,  very  accurate  results  could  be  obtained  at  comers  (GS,  OSG). 

The  initial  mathematical  singularity  of  the  solution  and  its  treatment  in  the  model  are 
discussed  in  the  next  Section.  Various  numerical  tests  are  then  carried  out  to  assess  accuracy  of 
results  obtained  using  new  corner  compatibility  conditions  (27)  versus  the  old  condition  (25), 
new  surface  interpolation  methods.  Mil  and  MCI,  versus  QS,  and  new  regridding  techniques. 
It  should  be  pointed  out  that,  in  the  following,  the  purpose  of  regridding  is  to  improve  result 
accuracy  in  two  situations  :  (i)  in  a  given  discretization,  by  adding  nodes  to  insufficiently 
resolved  free  surface  regions  (e.g.,  close  to  the  wavemaker);  and  (ii)  in  a  coarse  discretization,  by 
selectively  adding  nodes  to  specific  boundary  segments  (e.g.,  wavemaker,  wave  crest),  thereby 
obtaining  results  comparable  to  those  in  a  finer  discretization,  with  a  reduced  computational 
effort 14 . 


(i)  Wavemaker  motion  and  initial  boundary  conditions 

When  waves  are  generated  from  a  state  of  rest  using  a  wavemaker,  there  is  an  initial  jump  to 
finite  values  of  the  wavemaker  velocity  and  acceleration,  at  the  first  time  step  of  computations 
(on  boundary  Tr| .  Fig.  1),  whereas  initial  potential  c/>(x,  0)  and  elevation  rj (x,  0)  are  normally  set 
to  zero  on  the  free  surface  (boundary  T/).  This  creates  conflicting  boundary  conditions  on  both 
sides  of  a  comer  which,  in  fact,  are  similar  to  the  step  velocity  case  studied  by  Roberts  (1987) 
based  on  a  small  time  weakly  nonlinear  analysis.  Roberts’  results  showed  that  conflicting  initial 
boundary  conditions  create  a  singularity  associated  with  fast  small  scale  oscillations  of  the  free 
surface  close  to  the  wavemaker  15 .  The  same  study  also  predicted  non  singular  behavior  for  an 

14This  situation  is  similar  to  the  shoaling  computations  in  Section  3.2 

15Note  that,  due  to  these  fast  oscillations,  nonlinear  and  surface  tension  effects  should  play  a  role  and  should 
also  be  included  in  the  analysis.  This  was  somewhat  improved  in  the  solution  of  the  same  problem  by  Joo,  Schultz 
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initial  wavemaker  velocity  starting  from  zero  and  exponentially  increasing  with  time.  The  latter 
result  should  be  expected  since,  in  this  case,  both  the  initial  velocity  and  the  acceleration  are 
zero  and  the  initial  boundary  conditions  are  well  posed.  The  disadvantage  of  such  an  approach, 
however,  as  pointed  out  in  OSG,  is  that  it  would  take  an  infinite  time  for  the  wavemaker  to  reach 
a  finite  velocity. 


With  these  observations  in  mind,  GS  studied  the  accuracy  of  the  solution  close  to  the 
wavemaker,  when  generating  first-order  (Boussinesq)  solitary  waves  with  a  piston  wavemaker 
moving  as, 

H  k  k\ 

xP{t )  =  — -[tanh  — (ct  —  xp(t)  —  A)  +  tanh  — —  ]  (28) 

K/  /I'q  tip 

in  which  Ha  is  the  wave  height  in  constant  depth  hQ,  k  =  and  the  celerity  c  = 

\J g(ha  +  H0).  Piston  motion  (28)  corresponds  to  a  wave  profile  truncated  at  distance  A  on  both 
sides  of  the  wave  crest,  where  free  surface  amplitude  has  reduced  to  ezHa.  Developing  (28), 
the  initial  wavemaker  acceleration  reads, 


u. 


,(0)  =  V3gez(H’f2(l  +  H’o) 


1  ~  £z 


[l  +  ezHt 
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(29) 


which,  for  small  zz  and  fixed  H'a,  is  roughly  proportional  to  ez  and  can  thus  be  made  arbitrarily 
small  (initial  velocity  is  also  roughly  proportional  to  ez) 16 .  Thus,  for  small  ez  (i.e.,  for  large  A), 
both  the  initial  velocity  and  acceleration  are  small  but  finite  and  the  piston  motion  subsequently 
follows  a  smooth  “tanh”  law  (with  initial  exponential  increase). 


When  generating  solitary  waves  using  a  small  ez  value  (typically  0.002)  corresponding 
to  an  initial  acceleration  much  smaller  than  gravity,  GS  obtained  fairly  smooth  free  surface 
profiles  close  to  the  wavemaker  and  concluded  that,  for  all  practical  purposes,  limiting  the 
initial  acceleration  to  a  small  fraction  of  g  was  sufficient  to  eliminate  initial  singularity  problems 
l7.  In  a  more  detailed  analysis  of  these  results,  however,  OSG  showed  that  small  oscillations 
do  occur  on  the  free  surface  close  to  the  wavemaker,  first  in  the  higher-order  derivatives  of 

a2i  q2  j. 

the  potential  (g^,  and  propagate  downstream  with  the  generated  wave.  This  is  a  clear 
indication  of  initial  singularity  problems.  In  the  model,  however,  it  takes  a  fairly  large  distance 
of  propagation  for  these  oscillations  to  grow  and  to  significantly  affect  the  wave  profile.  In  the 
light  of  theoretical  conclusions  above,  OSG  re-calculated  the  generation  of  a  fairly  large  wave 
(H'0  =  0.4),  using  a  much  smaller  ez  =  3.5  10-5  (A  =  5.8 h0/ k)  corresponding  to  a  very  small 
initial  acceleration,  •Up(O)  =  2.2  10_5g.  Doing  so,  initial  inconsistencies  in  boundary  conditions 
were  very  much  reduced  and  results  indeed  showed  that  oscillations  in  higher-order  derivatives 
were  eliminated.  Initiation  of  the  wavemaker  motion  and  wave  generation,  however,  took  a 
much  larger  time  than  when  using  a  larger  ez. 

and  Messiter  (1990)  which  included  surface  tension  effects. 

l6For  small  ez ,  iip( 0)  ~  V3ez  (H'0)3/2(  1  +  H'0)  g  (with,  e.g.,  A  =  3.8/i0/re  for  ez  =  0.002). 

17This  also  qualitatively  agreed  with  Roberts’  (1987)  conclusions  when  the  exponent  in  his  power  law  was 
greater  than  2. 
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To  be  able  to  use  larger  ez  values  and  thus  reduce  the  wave  generation  time,  OSG  proposed  to 
initiate  their  calculations  with  consistent  boundary  conditions  on  the  free  surface,  by  specifying 
both  a  small  initial  elevation  and  potential  on  the  free  surface  corresponding  to  the  initial  velocity 
itp(O)  and  acceleration  Ap(0)  of  the  wavemaker.  To  do  so,  they  assumed  exponential  forms  for 
rj(x,  0)  and  (f>(x,  0)  and,  for  small  initial  amplitudes,  they  solved  the  corresponding  linearized 
problem  and  obtained, 

r,(x,0)=^e-*~  ,  ^,0)  =  -^e-2“  (30) 

Using  (30)  as  initial  conditions  on  the  free  surface  and  re-calculating  cases  above,  they  showed 
that  oscillations  could  be  eliminated  for  H'a  =  0.4,  even  with  ez  =  0.002,  which  confirmed  that 
initial  (mathematical)  singularities  were  effectively  removed.  This  method  is  used  in  the  present 
computations. 

When  eliminating  the  initial  (mathematical)  singularity  as  described  above,  the  impulsive 
flow  induced  by  a  strongly  moving  boundary  constitutes  a  very  demanding  numerical  test  case 
for  the  accuracy  of  compatibility  conditions  in  the  free  surface  corner  and  also  one  that  requires 
the  discretization  in  this  area  to  be  sufficiently  fine,  due  to  large  flow  variations  close  to  the 
comer,  and  the  BEM  integrals,  particularly  the  quasi-singular  ones,  to  be  calculated  with  great 
care  and  accuracy.  In  the  study  by  GS,  a  very  large  wave  height,  H'a  =  2,  was  selected  in  (28) 
(with  ez  =  0.002)  in  order  to  generate  a  wave  that  rapidly  overturns  and  breaks  quite  close  to 
the  wavemaker  18.  This  case  is  re-computed  here  in  order  to  test  the  accuracy  of  both  the  new 
interpolation  methods  and  new  comer  compatibility  condition,  versus  older  methods  used  in 
GS.  Eqs.  (30)  are  used  as  discussed  above  to  eliminate  initial  singularities  in  the  computations 
19.  After  carrying  out  these  calculations  in  the  same  two  discretizations  as  in  GS,  various  node 
regridding  strategies  are  tested  in  order  to  improve  both  the  accuracy  and  the  stability  of  results 
in  the  coarser  discretization. 


(ii)  Numerical  data  and  preliminary  computations 

The  computational  domain  is  rectangular  (as  in  Fig.  2)  with  a  depth  hD  =  1.  Two  spatial 
discretizations  are  used  in  the  computations  :  (i)  a  coarse  one  with  node  spacing  AsJ,  =  0.25 
on  the  free  surface  and  on  the  bottom  and  a  domain  length  10 hQ\  and  (ii)  a  fine  one  with 
AxJ,  =  0. 10  on  the  free  surface  and  on  the  bottom,  and  a  domain  length  8 hQ  2<).  In  the  first  case, 

18In  this  case,  although  maximum  wavemaker  acceleration  reached  one  g,  the  initial  acceleration  given  by  (29) 
was  still  small,  only  about  0.03 g. 

19Note,  however,  in  the  present  case,  due  to  the  short  propagation  distance  (and  propagation  time)  of  generated 
waves,  initial  small  oscillations  in  <ptn-  <ASS,  obtained  for  ez  =  0.002  when  initialization  (30)  is  not  used,  would 
not  significantly  affect  global  accuracy  of  computations. 

20The  shorter  computational  domain  was  selected  in  the  finer  discretization  to  reduce  computational  time.  This 
will  not  influence  results  at  all  in  the  region  up  to  and  including  the  breaker  jet,  in  the  last  computed  profile. 
Reflection  will  only  slightly  increase  in  the  wave  front,  for  x '  >  7  or  so,  and  will  lead  to  small  discrepancies  in 
free  surface  elevation  between  results  computed  in  both  discretizations. 
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A z'Q  =  0.125  on  lateral  boundaries  and,  in  the  second  case,  A z'Q  =  0.10.  The  total  number  of 
nodes  is  Nr  =  100  and  184,  for  each  case  respectively.  As  in  GS,  two-node  linear  elements 
are  used  on  lateral  boundaries  and  on  the  bottom,  and  adaptive  integration  is  specified  in  all 
comer  elements.  With  these  data,  the  computation  time  is  0.80  and  1.90  sec  per  time  step  (IBM 
9000/3),  for  each  discretization  respectively.  The  free  surface  interpolation  and  the  size  of  the 
time  step  are  discussed  in  the  following. 

A  large  wave  is  generated  in  both  discretizations  by  moving  the  lateral  boundary  according 
to  (28),  with  H'a  =  2.  Fig.  6  shows  free  surface  profiles  computed  with  the  MCI  method 
(with  C0  =  0.40),  using  the  new  compatibility  condition  (27)  in  free  surface  corners.  As 
expected,  due  to  the  strong  movement  of  the  lateral  boundary,  the  generated  wave  overturns 
with  a  large  scale  plunging  jet  21 .  Figs.  6a-b  show  results  in  the  finer  discretization  and  Fig. 
6c  shows  results  in  the  coarser  discretization.  Fig.  6d  is  a  comparison  between  results  in 
both  discretizations.  In  the  coarser  discretization,  likely  due  to  insufficient  resolution  in  the 
breaker  jet,  computations  become  inaccurate  and  have  to  be  stopped  at  an  earlier  stage  (curve 
i)  than  in  the  finer  discretization  (curve  j)  where  computations  can  accurately  be  pursued  up  to 
impact  of  the  breaker  jet  on  the  free  surface  (the  theoretical  limitation  of  FNPF  models).  The 
latter  is  a  marked  improvement  versus  computations  in  GS,  using  the  QS  method  in  the  finer 
discretization,  which  had  to  be  stopped  at  the  stage  of  curve  f.  Accuracy  of  computations  in  the 
coarser  discretization  will  be  improved  by  regridding  (Section  v  below).  Fig.  6d  shows  that,  for 
curves  a-i,  results  obtained  with  both  discretizations  are  in  quite  good  agreement,  except  close 
to  the  wavemaker  and  in  the  jet  regions  (assuming  reflection  in  the  wave  front  is  ignored). 


(iii)  Comparison  between  the  old  and  the  new  corner  compatibility  conditions 

To  test  the  old  comer  compatibility  condition  (25)  used  in  GS  (subscript  o)  versus  new  condition 
(27)  (subscript  n),  computations  are  carried  out  in  both  discretizations  using  the  QS  interpolation 
method  on  the  free  surface,  with  Ca  =  0.32  (subscript  1)  and  0.40  (subscript  3).  Relative  errors 
in  wave  mass22  are  given  in  Fig.  7a  as  curves  (QSi0,  QSin)  and  (QS3o,  QS3n).  Comparing 
these  curves  two-by-two  for  both  discretizations,  it  is  clear  that  using  new  comer  compatibility 
conditions  leads  to  more  accurate  results  and  to  longer  propagation  times  before  instability 
occurs  than  using  the  old  conditions. 


(iv)  Comparison  between  the  QS  and  the  Mil  and  MCI  interpolation  methods 

To  assess  the  accuracy  of  new  interpolation  methods  on  the  free  surface,  computations  are 
carried  out  using  the  QS,  Mil,  and  MCI  interpolation  methods,  with  the  new  comer  compatibility 
condition  (27)  and  a  mesh  Courant  number  equal  to  0.40  (subscript  3). 

21Note  the  bimodal  tip  of  the  breaker  jet,  resulting  from  the  somewhat  unusual  way  of  generating  breaking  with 
a  piston  wavemaker. 

22with  m  =  3.26,  the  maximum  wave  mass  above  the  undisturbed  free  surface  z  =  0. 
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Relative  errors  in  wave  mass  are  given  in  Fig.  7a  for  these  computations.  They  show 
that  results  obtained  with  the  Mil  method,  although  initially  almost  identical  to  those  obtained 
with  the  MCI  method,  quickly  become  unstable  in  both  discretizations  whereas  results  with 
the  MCI  method  can  accurately  be  calculated  for  a  longer  time.  In  the  finer  discretization,  the 
instability  with  the  Mil  method  occurs  shortly  after  the  wave  reaches  the  breaking  point  (curve 
e  in  Fig.  6),  hence,  even  before  the  breaker  jet  starts  forming,  indicating  that  the  instability  is 
not  due  to  problems  in  the  breaker  jet.  More  likely,  as  nodes  move  away  from  the  wavemaker 
(due  to  node  drift),  results  in  the  first  interval  of  the  first  Mil  element  on  the  free  surface  are 
getting  increasingly  poor.  With  time,  this  leads  to  an  accumulation  of  errors  in  the  comer  and, 
eventually,  to  the  instability  of  computations.  The  spline-based  MCI  method,  on  the  other  hand, 
does  not  extrapolate  the  geometry  close  to  the  comer  and,  hence,  is  not  so  sensitive  to  a  loss 
of  resolution  of  the  discretization.  As  expected  from  shoaling  computations  in  Section  3.2, 
MCI  results  are  also  more  accurate  than  those  obtained  with  the  QS  method.  These  conclusions 
justify  the  implementation  of  the  MCI  method  for  problems  with  moving  solid  boundaries. 

The  influence  of  the  mesh  Courant  number  on  the  accuracy  of  results  with  the  MCI  method 
is  illustrated  in  Fig.  7b.  In  both  discretizations,  the  smallest  errors  and  the  longest  propagation 
times  are  obtained  using  the  largest  Courant  number,  CQ  =  0.40.  In  fact,  as  discussed  in  Section 
3.1,  using  too  small  a  Courant  number  leads  to  an  unnecessary  large  numbers  of  time  steps  and 
to  an  accumulation  of  numerical  errors. 

Finally,  whereas  relative  errors  on  wave  volume  at  the  end  of  calculations  are  on  the  same 
order  in  both  discretizations  (~0.3%,  in  Figs.  7a-b),  energy  errors  are  more  than  three  times 
larger  in  the  coarser  discretization  than  in  the  finer  discretization  (Fig.  7c).  This  is  because  the 
poorer  resolution  on  the  free  surface  in  the  coarser  discretization  leads  to  less  accurate  results 
for  the  kinematics  and,  hence,  for  the  wave  energy. 


(v)  Node  regridding  close  to  the  wavemaker  and  in  the  breaker  jet 

In  previous  computations,  numerical  errors  mainly  resulted  from  a  lack  of  resolution  of  the 
discretization  close  to  the  wavemaker,  due  to  node  drift  (Fig.  6b-c,  curves  a-j),  and  from 
quasi-singular  situations  created  in  breaker  jets,  due  to  node  proximity.  Regridding  will  now 
be  used  to  selectively  add  nodes  close  to  the  wavemaker,  on  the  free  surface,  and  adaptive 
regridding  will  be  used  to  prevent  nodes  from  moving  too  close  to  each  other  in  breaker  jets. 
Various  regridding  strategies  addressing  these  problems  (Table  1,  a-i)  were  tested  in  the  coarser 
discretization  (AsJ,  =  0.25),  using  the  MCI3  method,  and  corresponding  numerical  errors  and 
wave  profiles  are  given  in  Figs.  8  and  9,  respectively.  Time  t\  =  3.93  for  the  first  regridding 
corresponds  to  the  breaking  point  (vertical  wave  front,  i.e.,  approximately  curve  e  in  Fig.  6). 
Errors  in  Fig.  8  give  a  measure  of  the  effectiveness  of  each  regridding  method  in  allowing 
computations  to  be  carried  out  in  the  coarser  discretization,  with  sufficient  accuracy,  for  a  time 
longer  than  without  regridding  (curve  MCI3). 
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In  Fig.  8a,  strategies  a  or  b  (i.e.,  a  simple  regridding  of  the  whole  free  surface  without/with 
adjustment  of  CQ)  lead  to  increased  errors  without  allowing  computations  to  be  carried  out  for 
a  longer  time.  This  is  likely  because  breaker  jets  are  regridded  to  too  large  a  node  interval. 
Strategy  c  (i.e.,  regridding  of  just  the  first  7  nodes)  does  not  initially  increase  errors — since  the 
breaker  jet  is  not  regridded — but  leads  to  instability  of  computations  at  an  early  time.  Strategy 
d  (i.e.,  the  addition  of  10  nodes  close  to  the  wavemaker  at  a  later  time  t'2  after  regridding  of 
the  whole  free  surface  at  t\)  allows  computations  to  be  carried  out  for  a  slightly  longer  time 
than  without  regridding,  but  final  errors  are  more  than  three  times  as  large.  Strategy  e  (i.e., 
the  addition  of  10  nodes  close  to  the  wavemaker  at  time  t\)  leads  to  more  accurate  results  for 
most  of  the  propagation  (the  error  curve  in  Fig.  8a  is  almost  flat),  but  computations  cannot  be 
carried  out  for  a  longer  time  than  without  regridding.  Fig.  9a  shows  that  the  last  computed 
wave  profile  in  this  case  (curve  B)  is  in  good  agreement  with  results  in  the  finer  discretization. 
Some  discrepancies  are  still  observed  close  to  the  wavemaker  but  are  smaller  than  without 
regridding  (see  Fig.  6d).  Strategy  f,  (a  variation  on  strategy  e  that  adds  10  more  nodes  close  to 
the  wavemaker  at  a  later  time  t'2)  does  not  further  improve  the  accuracy. 

In  Fig.  8b,  strategies  g  or  h  (performing  additional  regridding  of  the  first  30  nodes  including 
the  breaker  jet  at  time  t'2  with/without  adjustment  of  Ca)  allow  computations  to  be  carried  out 
for  a  time  almost  as  long  as  in  the  finer  discretization  it'  =  5.35)  but  with  slightly  larger  errors. 
Figs.  9b-c  show  that  the  last  computed  wave  profiles  in  these  cases  (curves  B)  are  in  good 
agreement  with  results  in  the  finer  discretization.  In  strategies  g  and  h,  regridding  of  the  breaker 
jet  at  time  t'2  actually  reduces  the  tendency  of  nodes  to  move  too  close  to  each  other  and  thus 
also  reduces  quasi-singular  situations.  A  better  way  of  achieving  this  is  to  use  the  adaptive 
regridding  method  that  checks  node  distances  at  all  time  steps  and  regrids  nodes  two-by-two  if 
needed.  This  is  finally  done  in  strategy  i,  which  combines  strategy  e  to  improve  resolution  close 
to  the  wavemaker  and  the  adaptive  regridding  to  improve  accuracy  in  the  breaker  jet.  Error 
curves  in  Fig.  8  for  strategy  i  clearly  give  the  best  accuracy  and  the  longest  propagation  time 
of  all  regridding  methods — in  fact  as  long  as  in  the  fine  discretization — .  Final  error  on  wave 
mass  is  similar  to  the  error  in  the  finer  discretization  but  the  error  on  total  energy  is  larger.  Fig. 
9d  shows  that  wave  profiles  computed  with  strategy  i  agree  quite  well  with  those  obtained  in 
the  finer  discretization,  up  to  impact  of  the  jet  (curve  j).  Only  the  smaller  scale  features,  like 
the  bimodal  breaker  jet,  are  not  represented  due  to  insufficient  resolution.  Comparing  curves  e 
in  Fig.  6c  and  9d,  one  also  sees  how  regridding  has  added  nodes  in  the  back  of  the  wave  and, 
comparing  curves  f-i  in  both  figures,  how  nodes  have  now  moved  much  less  apart,  close  to  the 
wavemaker. 

Finally,  it  should  be  stressed  that  results  obtained  with  strategy  i  are  calculated  in  a  dis¬ 
cretization  that  has  only  10  more  nodes  than  the  original  coarse  discretization  (total  110).  This 
leads  to  a  slight  increase  in  computation  time  to  about  0.95  sec  per  time  step,  i.e.,  only  half  the 
time  needed  in  the  finer  discretization.  Hence,  using  the  regridding  techniques,  accurate  results 
can  be  obtained  in  a  coarser  discretization  with  a  much  reduced  computational  effort. 
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4  Conclusions 


New  interpolation  methods  and  corner  compatibility  conditions  were  implemented  in  an  existing 
fully  nonlinear  wave  propagation  model  based  on  a  BEM.  Results  showed  improvements  in  both 
accuracy  and  stability  of  computations  for  wave  breaking  induced  by  fixed  or  moving  boundaries, 
as  compared  to  earlier  methods.  These  improvements  are  achieved  within  essentially  the  same 
CPU  time  (i.e.,  computational  effort)  as  with  the  earlier  model.  New  numerical  methods  are 
thus  also  found  to  be  more  efficient. 

Node  regridding  techniques  were  implemented  to  further  improve  both  the  accuracy  and 
the  duration  of  computations  before  numerical  instabilities  occur.  These  techniques  selectively 
add  and  redistribute  nodes  over  specified  sections  of  the  free  surface  and/or  adaptively  regrid 
nodes  two-by-two  on  the  free  surface  when  they  move  too  close  to  each  other.  The  combina¬ 
tion  of  regridding  techniques  with  new  interpolation  and  corner  compatibility  methods  allows 
computations  of  wave  breaking  (induced  by  a  bottom  slope  or  by  a  moving  boundary)  to  be 
accurately  pursued  beyond  the  breaking  point  (vertical  front  slope  of  the  wave),  up  to  impact  of 
the  breaker  jet  on  the  free  surface  (theoretical  limitation  of  FNPF  theory). 

With  the  earlier  model,  calculations  in  large  computational  domains,  in  particular,  had  to 
be  stopped  at  or  shortly  after  the  breaking  point  was  reached  thus  preventing  a  close  analysis 
of  post-breaking  wave  characteristics.  This  is  now  possible  using  the  new  methods  and  present 
results,  for  instance,  show  the  detailed  development  of  the  plunging  breaking  of  solitary  waves 
over  a  gentle  slope.  Based  on  such  computations,  detailed  internal  kinematics  of  breaking  waves 
over  a  slope  (or  over  an  arbitrary  bottom  topography)  can  be  calculated  (e.g.,  Grilli,  Svendsen, 
Subramanya  (1994);  Grilli  and  Horrillo  (1995)).  To  our  knowledge,  such  results  for  breaking 
waves  have  never  been  obtained  through  numerical  modeling. 
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Figure  captions 


Figure  1  :  Sketch  of  typical  computational  domain  for  solitary  wave  shoaling  over  a  gentle  1:35 
slope. 


Figure  2  :  Sketch  of  computational  domain  for  the  propagation  of  a  fully  nonlinear  solitary 
waves  of  height  H/h  =  0.3  —  0.6,  over  5  units  of  time  (domain  length  is  L/h  =  28).  Results 
are  given  in  Fig.  3. 


Figure  3  :  Absolute  numerical  errors  (A’s)  with  respect  to  theoretical  values  by  Tanaka  :  (abc) 
for  H/h  =  0.5;  or  (d)  as  a  function  of  H/h.  In  (abc),  errors  in  mass  m,  energy  e,  and  height  H, 

are  given  as  a  function  of  the  free  surface  interpolation  method  :  ( - )  MCI;  ( - )  Mil;  and 

( - )  QS,  and  the  initial  distance  between  nodes  :  Ax'a  =  (□)  0.25;  (o)  0.20;  and  (o)  0.15.  In 

(d),  errors  are  given  for  Ax'a  =  0.15  (MCI  interpolation)  and  CQ  =  0.35  (adaptive  time  step). 


Figure  4  :  Blow-ups  of  Fig.  1  for  the  shoaling  and  breaking  of  solitary  waves  of  initial  height 
:  (ab)  H'0  =  0.20;  and  (cd)  0.40,  over  a  1:35  slope.  Curves  :  a-d  and  A-D,  are  without  regrid- 
ding;  e-k  and  E-I,  are  after  regridding.  (o)  denote  discretization  nodes,  (-  -  o  -  -)  represent  last 
computed  profile  without  regridding,  and  ( — o — )  represent  regridded  profiles.  Time  of  plotted 
wave  profiles  is  t'  =  a  :  37.17,  b  :  40.73,  c  :  43.48,  d  :  44.53,  e  :  44.62,  f  :  44.94,  g  :  45.20, 
h  :  45.40,  i  :  45.60,  j  :  45.80,  k  :  46.00,  A  :  13.15,  B  :  16.00,  C  :  17.67,  D  :  18.33,  E  :  18.66, 
F :  19.00,  G :  19.30,  H :  19.60,  andl :  19.90.  Corresponding  numerical  errors  are  plotted  in  Fig.  5. 


Figure  5  :  Shoaling  of  solitary  waves  of  initial  height  H'a  over  a  1:35  slope  :  (a),(d)  relative 
errors  on  wave  volume  m;  (b)  relative  errors  on  total  energy  e;  and  (c)  shoaling  curves  H/h. 

Free  surface  interpolation  method  :  ( - )  MCI;  ( - )  QS;  ( - )  MCI  with  regridding. 

Arrows  in  Fig.  5a  indicate  locations  at  which  errors  on  total  energy  reach  1.0%  in  Fig.  5b. 


Figure  6  :  Wave  breaking  induced  by  a  piston  wavemaker  (solitary  wave  motion  with  if'  =  2.0) 
with  :  Ax'0  and  CD  =  (a)(b)  0.10,  0.40;  (c)  0.25,  0.40.  (d)  shows  superimposed  results  of  : 

( - )  (a),(b)  and  ( - )  (c).  (o)  denote  discretization  nodes  and  vertical  lines  mark  successive 

locations  of  the  piston  wavemaker.  Time  of  plotted  wave  profiles  is  t'  =  a  :  1.442,  b  :  2.064,  c 
:  2.843,  d:  3.381,  e:  4.011,  f:  4.353,  g  :  4.681,  h:  4.936,  i :  5.191,  j:  5.438. 


Figure  7  :  Relative  errors  in  mass  m  and  total  energy  e  for  the  calculations  in  Fig.  6  with  Ax  'Q  = 

( - )  0.10;  ( - )  0.25,  for  three  interpolation  methods:  QS,  Mil,  and  MCI.  Subscripts  refer 

to  :  (i)  C0  =  (1)  0.32;  (2)  0.35;  (3)  0.40;  (ii)  Node  compatibility  =  (o)  old  method;  (n)  new 
method  (also  used  if  nothing  else  is  mentioned). 
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Figure  8  :  Relative  errors  in  the  coarser  discretization,  for  the  same  calculations  as  in  Fig. 
6-7,  using  regridding  strategies  (a-i)  in  Table  1  with  MCI3  and  AxJ,  =  0.25.  Curves  MCI3 
correspond  to  results  obtained  without  regridding. 


Figure  9  :  Free  surface  profiles  obtained  for  the  same  calculations  as  in  Fig.  8  using  regridding 

method  :  (a)  e;  (b)  g;  (c)  h;  and  (d)  i,  compared  to  computations  in  the  fine  discretization  ( - ) 

with  MCI3  (no  regridding,  also  shown  in  Fig.  6a).  (o)  denote  discretization  nodes  and  results. 
Time  of  plotted  wave  profiles  is  t'  =  (a)  A  :  5.05,  B  :  5.25;  (b)  A  :  5.05,  B  :  5.30;  (c)  A  :  5.05,  B 
:  5.35;  (d)  e-j  as  in  Fig.  6. 
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Strategy 

Description  of  regridding  strategy 

a 

At  t  =  t\,  regridding  of  the  free  surface  to  constant  step. 

b 

At  t  =  1 1 ,  regridding  of  the  free  surface  to  constant  step,  CQ  =  0.40. 

c 

Att  =  t\,  regridding  of  the  first  7  nodes  to  constant  step. 

d 

Same  as  b  and,  at  t  =  t2,  add  and  regrid  10  nodes  to  the  first  6  node  segment. 

e 

At  t  =  fi,  add  and  regrid  10  nodes  to  the  first  6  node  segment. 

f 

Same  as  e  and,  at  t  =  ti,  add  and  regrid  10  nodes  to  the  first  6  node  segment. 

g 

Same  as  e  and,  at  t  =  t2,  regrid  first  30  node  segment  to  constant  step,  CQ  =  0.40. 

h 

Same  as  e  and,  at  t  =  t2,  regrid  first  30  node  segment  to  constant  step. 

i 

Same  as  e  and,  for  t  >  t\ ,  adaptive  regridding  on  the  free  surface. 

Table  1:  Regridding  strategies  for  the  computations  in  Fig.  8  and  9,  with  MCI3  and  Ax'a  =  0.25 
(coarse  grid).  Regridding  times  :  t\  =  3.93,  and  t'2  =  5.00. 
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